Initialisation
## [1] "number of cores available = 16"
#Phi[1] ; eta = valeur de fin
#Phi[2] = valeur du noeud
#Phi[3] = echelle
m <- function(t, eta, phi) (phi[1] + eta)/(1+exp((phi[2]-t)/phi[3]))
#=======================================#
param <- list(sigma2 = 0.05,
rho2 = 0.1,
mu = c(5,90,5),
omega2 = c(0.5,0.1,0.01),
#Survival data,
nu2 = 0.5,
a = 90,
b = 50,
alpha = 7,
beta = 10)
F. <- length(param$omega2) #dimension de phi
#=======================================#
t <- seq(60,120, length.out = 10) #value of times
# --- longitudinal data --- #
G = 40 ; ng = 4 #nombre de groupe et d'individu par groupe
dt_NLME <- NLME_obs(G, ng, t, param, m)
Y <- dt_NLME$obs
dt_NLME %>% ggplot(aes(t, obs, col = gen, group = id)) +
geom_point() + geom_line() + theme(legend.position = 'null')

# --- Survival data --- #
dt_SF <- SF_obs(dt_NLME, param, m)
dt_SF %>% ggplot(aes(obs, fill = factor(gen))) +
geom_histogram(position = 'dodge', bins = 20) + theme(legend.position = 'null')

SAEM avec simulation par MCMC
\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]
sigma2_a <- sigma2_b <- sigma2_alpha <- 0.1
Phi <- function(sigma2, rho2, mu, omega2, nu2, a, b, alpha, beta)
c(- 1/(2*sigma2), -1/(2*rho2), -1/(2*omega2), G*mu/omega2, #1, 2, 3, 4
-1/(2*nu2), a/sigma2_a, -1/(2*sigma2_a), #5, 6, 7
b/sigma2_b , -1/(2*sigma2_b), #8, 9
alpha/sigma2_alpha, -1/(2*sigma2_alpha) ) #10, 11
zeta <- function(beta, eta, phi, gamma, a, b, alpha)
{
lbd <- function(t,g) t^{b-1} * exp(alpha*m(t, eta[g], phi[g,]))
P1 <- 1:nrow(dt_SF) %>% sapply(function(i) integrate(function(t) lbd(t,dt_SF$gen[i]) , 0, dt_SF$obs[i])$value )
P2 <- 1:nrow(dt_SF) %>% sapply(function(i) dt_SF$U[i]*exp(beta*dt_SF$U[i] + gamma[dt_SF$gen[i]] ) )
sum(dt_SF$U) - b*a^-b * sum(P2*P1)
}
S <- function(eta, phi, gamma, a, b, alpha)
{
beta <- uniroot(function(beta)zeta(beta, eta, phi, gamma, a,b,alpha), lower = -100, upper = 100)$root
s <- c(mean((Y - get_obs(m, dt_NLME, eta = eta, phi = phi) )^2 ), #1
mean(eta^2), #2
apply(phi^2, 2, mean), apply(phi , 2, mean), #3 et 4
mean(gamma^2), #5
a, a^2, #6 et 7
b, b^2, #8 et 9
alpha^2, alpha, #10 et 11
beta ) #12
attr(s, '1') <- 1
attr(s, '2') <- 2
attr(s, '3') <- 2 + 1:ncol(phi)
attr(s, '4') <- 2 + ncol(phi) + 1:ncol(phi)
attr(s, '5') <- 2 + 2*ncol(phi) + 1
for(i in 6:12) attr(s, as.character(i)) <- attr(s, as.character(i-1)) + 1
return(s)
}
#=============================================================================#
loglik.phi <- function(x, eta, Phi)
{
Phi[1] * sum((Y - get_obs(m, dt_NLME, eta = eta, phi = x) )^2) +
Phi[3] * apply(x^2, 2, sum) + Phi[4] * apply(x , 2, sum)
}
loglik.eta <- function(x, phi, Phi)
{
Phi[1] * sum((Y - get_obs(m, dt_NLME, eta = x, phi = phi) )^2) +
Phi[2] * sum(x^2)
}
loglik.gamma <- function(x, Phi) Phi[5] * sum(x^2) + ng*sum(x) #Correction fait à la va vite
loglik.a <- function(x, b, Phi) Phi[6] * x + Phi[7] * x^2 + G*ng*log(b*x^-b)
loglik.b <- function(x, Phi) Phi[8] * x + Phi[9] * x^2 + (x-1)*sum(log(t))
loglik.alpha <- function(x, Phi) Phi[10] * x + Phi[11] * x^2
SAEM
initialisation
M <- 1 #nombre de simulation
u <- function(k) ifelse(k<200, 1, 1/(k-199) )
# --- Initialisation des paramètres --- #
para <- param %>% lapply(function(x) x + 2*runif(1) )
# --- Initialisation des chaines MC : Z_0 --- #list(attr(dt_NLME,'eta')),#
Z <- list(eta = 1:M %>% lapply(function(i) rnorm(G*ng, 0, para$rho2) %>% matrix(ncol = 1) ),
phi = 1:M %>% lapply(function(i) matrix(rnorm(F.*G, para$mu, para$omega2), nrow = F.) %>% t ),
#list(attr(dt_NLME,'phi')),#
gamma = 1:M %>% lapply(function(i) matrix(rnorm(G, 0, para$nu2), ncol = 1) ),
a = list(matrix(para$a)), #list(matrix(param$a)), #
b = list(matrix(para$b)), #list(matrix(param$b)), #
alpha = list(matrix(para$alpha)) )#list(matrix(param$alpha)) )#
Boucle
sim <- function(niter, Phih, eta, phi, gamma, a, b, alpha)
{
M <- length(phi)
eta <- 1:M %>% lapply( function(i)
MH_High_Dim_para(niter, eta[[i]], sd = .02, loglik.eta, phi[[i]], Phih, cores = ncores-1 ))
# | | |
phi <- 1:M %>% lapply( function(i) # | |
MH_High_Dim_para(niter, phi[[i]], sd = c(.03, .02, .02), loglik.phi, eta[[i]], Phih, cores = ncores-1 ))
# | | |
gamma <- 1:M %>% lapply( function(i)# | |
MH_High_Dim_para(niter, gamma[[i]], sd = .03, loglik.gamma, Phih, cores = ncores-1 ))
a <- 1:M %>% lapply( function(i) MH_High_Dim_para(niter, a[[i]], sd = .02, loglik.a,b[[i]], Phih, cores = 1 ))
b <- 1:M %>% lapply( function(i) MH_High_Dim_para(niter, b[[i]], sd = .02, loglik.b, Phih, cores = 1 ))
alpha <- 1:M %>% lapply( function(i) MH_High_Dim_para(niter, alpha[[i]], sd = .02, loglik.alpha, Phih, cores = 1 ))
list(eta = eta, phi = phi, gamma = gamma, a = a, b = b, alpha = alpha)
}
maxi <- function(S)
{
list(sigma2 = S[attr(S,'1')],
rho2 = S[attr(S,'2')],
mu = S[attr(S,'4')],
omega2 = S[attr(S,'3')] - S[attr(S,'4')]^2,
nu2 = S[attr(S,'5')],
a = S[attr(S,'6')], b = S[attr(S,'8')],
alpha = S[attr(S,'11')], beta = S[attr(S,'12')] )
}
Estimation oracle
|
sigma2
|
rho2
|
mu1
|
mu2
|
mu3
|
omega21
|
omega22
|
omega23
|
nu2
|
a
|
b
|
alpha
|
beta
|
|
0.0511813
|
0.1045811
|
4.932414
|
90.07442
|
4.980973
|
0.4093858
|
0.0762974
|
0.0081025
|
0.5462794
|
90
|
50
|
7
|
9.967896
|
niter <- 300
MH.iter <- 10
res <- SAEM(niter, MH.iter, u, para, Phi, S, Z, sim, maxi, eps = 1e-1, verbatim = T)
saveRDS(res, 'saem.rds')
gg <- plot(res)
## [1] "SAEM execution time = 00min 02sec"
Résultat de l’algo SAEM-MCMC
|
|
sigma2
|
rho2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
nu2
|
a
|
b
|
alpha
|
beta
|
|
Valeur réelle
|
0.0500
|
0.1000
|
5.0000
|
90.0000
|
5.0000
|
0.5000
|
0.1000
|
0.0100
|
0.5000
|
90.0000
|
50.0000
|
7.0000
|
10.0000
|
|
Valeur estimée
|
0.0517
|
3.4131
|
3.6631
|
89.9177
|
5.0576
|
1.7349
|
0.3362
|
0.2266
|
0.8383
|
107.1835
|
68.2192
|
23.7318
|
20.0341
|
|
Rrmse
|
0.0342
|
33.1310
|
0.2674
|
0.0009
|
0.0115
|
2.4697
|
2.3615
|
21.6611
|
0.6766
|
0.1909
|
0.3644
|
2.3903
|
1.0034
|


